Automobile Insurance Models

Introducción

Lecturas Sugeridas:

Supervised Machine Learning: Classification and Regression https://medium.com/@nimrashahzadisa064/supervised-machine-learning-classification-and- regression-c145129225f8

What is Supervised Learning? https://www.ibm.com/topics/supervised-learning

n A Beginner’s Guide to Supervised Machine Learning Algorithms https://towardsdatascience.com/a-beginners-guide-to-supervised-machine-learning-algorithms- 6e7cd9f177d5

Teoría de Modelos y Algoritmos

Brevemente responder con tus propias palabras 2 de las siguientes 3 preguntas: - i) ¿Qué es Supervised Machine Learning y cuáles son algunas de sus aplicaciones en Inteligencia de Negocios?

• 2. ¿Cuáles son los principales algoritmos de Supervised Machine Learning? Brevemente describir con tus propias palarbas 5 – 7 de los principales algoritmos de Supervised Machine Learning.
• 3. ¿Qué es la R2 Ajustada? ¿Qué es la métrica RMSE? ¿Cuál es la diferencia entre la R2 Ajustada y la métrica RMSE?

Librerías Necesarias

# Manipulación y Visualización de Datos
library(ggplot2)      # Sistema para crear gráficos
library(DataExplorer) # Facilita la exploración de datos
library(naniar)       # Facilita la visualización de valores faltantes y la examinación de imputaciones
library(dlookr)       # Herramientas para diagnóstico de datos, exploración y transformación
library(RColorBrewer) # Ofrece varias paletas de colores
library(foreign)      # Lee datos almacenados por diversas plataformas

# Visualización y Manipulación de Datos Espaciales
library(sf)           # Una forma estandarizada de codificar datos vectoriales espaciales
library(mapview)      # Crea visualizaciones interactivas de datos espaciales rápidamente
library(tmap)         # Para dibujar mapas temáticos

# Modelado Predictivo
library(regclass)     # Herramientas básicas para visualizar, interpretar y construir modelos de regresión
library(mctest)       # Diagnóstico de multicolinealidad
library(lmtest)       # Pruebas para modelos de regresión lineal
library(spdep)        # Funciones para crear objetos de matriz de pesos espaciales
library(spData)       # Diversos conjuntos de datos espaciales
library(spatialreg)   # Funciones de estimación para modelos espaciales
library(caret)        # Funciones para simplificar el proceso de entrenamiento de modelos
library(e1071)        # Variedad de funciones para análisis
library(SparseM)      # Funcionalidad básica para álgebra lineal con matrices dispersas
library(Metrics)      # Implementación de métricas de evaluación comúnmente usadas en aprendizaje supervisado
library(tidyr)
library(readr)
library(randomForest) # Basado en un bosque de árboles usando entradas aleatorias
library(jtools)       # Herramientas para entender y compartir resultados de análisis de regresión
library(xgboost)      # Incluye solucionadores de modelos lineales eficientes y algoritmos de aprendizaje de árboles
library(DiagrammeR)   # Construye estructuras de grafos/redes
library(effects)      # Muestra gráfica y tabular de efectos, como interacciones, para modelos estadísticos
library(shinyjs)      # Funciones de JavaScript para aplicaciones Shiny
library(sp)           # Clases y métodos para datos espaciales
library(geoR)         # Análisis geostadístico
library(gstat)        # Modelado geostadístico, predicción y simulación
library(corrplot)
library(tigris)
library(car)
library(lmtest)
library(neuralnet)
library(MASS) 
library(dplyr)        # Herramienta rápida y consistente para trabajar con objetos similares a data frames

Obtención de Data

automoble_insurance_claims <- read_csv("automoble_insurance_claims.csv")

## Rows: 1000 Columns: 40
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (21): policy_bind_date, policy_state, policy_csl, insured_sex, insured_e...
## dbl (18): months_as_customer, age, policy_number, policy_deductable, policy_...
## lgl  (1): _c39
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(automoble_insurance_claims)

## # A tibble: 6 × 40
##   months_as_customer   age policy_number policy_bind_date policy_state
##                <dbl> <dbl>         <dbl> <chr>            <chr>       
## 1                328    48        521585 10/17/2014       OH          
## 2                228    42        342868 6/27/2006        IN          
## 3                134    29        687698 9/6/2000         OH          
## 4                256    41        227811 5/25/1990        IL          
## 5                228    44        367455 6/6/2014         IL          
## 6                256    39        104594 10/12/2006       OH          
## # ℹ 35 more variables: policy_csl <chr>, policy_deductable <dbl>,
## #   policy_annual_premium <dbl>, umbrella_limit <dbl>, insured_zip <dbl>,
## #   insured_sex <chr>, insured_education_level <chr>, insured_occupation <chr>,
## #   insured_hobbies <chr>, insured_relationship <chr>, `capital-gains` <dbl>,
## #   `capital-loss` <dbl>, incident_date <chr>, incident_type <chr>,
## #   collision_type <chr>, incident_severity <chr>, authorities_contacted <chr>,
## #   incident_state <chr>, incident_city <chr>, incident_location <chr>, …

Analisis Exploratorio de Datos

Desarrollar Análisis Exploratorio de los Datos (EDA) que incluye los siguientes elementos: ## a. Identificación de NA’s

# Cambia "?" por NA en todo el dataframe
df_aic <- data.frame(lapply(automoble_insurance_claims, function(x) {
  if (is.factor(x)) x <- as.character(x)
  x[x == "?"] <- NA
  if (is.character(x)) x <- factor(x)
  return(x)
}))

# Cuenta la cantidad de NA en cada columna
na_count <- sapply(df_aic, function(x) sum(is.na(x)))

# Imprime el conteo de NA de manera estética
cat("Conteo de NA por columna en 'automoble_insurance_claims':\n")

## Conteo de NA por columna en 'automoble_insurance_claims':

for (columna in names(na_count)) {
  cat(columna, ": ", na_count[columna], "\n")
}

## months_as_customer :  0 
## age :  0 
## policy_number :  0 
## policy_bind_date :  0 
## policy_state :  0 
## policy_csl :  0 
## policy_deductable :  0 
## policy_annual_premium :  0 
## umbrella_limit :  0 
## insured_zip :  0 
## insured_sex :  0 
## insured_education_level :  0 
## insured_occupation :  0 
## insured_hobbies :  0 
## insured_relationship :  0 
## capital.gains :  0 
## capital.loss :  0 
## incident_date :  0 
## incident_type :  0 
## collision_type :  178 
## incident_severity :  0 
## authorities_contacted :  0 
## incident_state :  0 
## incident_city :  0 
## incident_location :  0 
## incident_hour_of_the_day :  0 
## number_of_vehicles_involved :  0 
## property_damage :  360 
## bodily_injuries :  0 
## witnesses :  0 
## police_report_available :  343 
## total_claim_amount :  0 
## injury_claim :  0 
## property_claim :  0 
## vehicle_claim :  0 
## auto_make :  0 
## auto_model :  0 
## auto_year :  0 
## fraud_reported :  0 
## X_c39 :  1000

b. Reemplazo de NA’s

# Definir rangos para 'injury_claim'
df_aic$injury_claim_range <- cut(df_aic$injury_claim, breaks = quantile(df_aic$injury_claim, probs = seq(0, 1, by = 0.25), na.rm = TRUE), include.lowest = TRUE)

# Función para aplicar operaciones a cada grupo, ajustada para trabajar con rangos de 'injury_claim'
apply_operations_per_group <- function(df) {
  for (col in c("police_report_available", "property_damage", "collision_type")) {
    moda <- names(sort(table(df[[col]]), decreasing = TRUE))[1]
    df[[col]][is.na(df[[col]])] <- moda
  }

  return(df)
}

# Segmentar por 'injury_claim_range' y aplicar la función definida
df_aic <- df_aic %>%
  group_by(injury_claim_range) %>%
  group_modify(~apply_operations_per_group(.x))

# Opcional: Puedes remover la columna 'injury_claim_range' después de la imputación si ya no la necesitas
df_aic$injury_claim_range <- NULL
df_aic$X_c39 <- NULL

na_count <- sapply(df_aic, function(x) sum(is.na(x)))
cat("Conteo de NA por columna en 'automoble_insurance_claims':\n")

## Conteo de NA por columna en 'automoble_insurance_claims':

for (columna in names(na_count)) {
  cat(columna, ": ", na_count[columna], "\n")
}

## months_as_customer :  0 
## age :  0 
## policy_number :  0 
## policy_bind_date :  0 
## policy_state :  0 
## policy_csl :  0 
## policy_deductable :  0 
## policy_annual_premium :  0 
## umbrella_limit :  0 
## insured_zip :  0 
## insured_sex :  0 
## insured_education_level :  0 
## insured_occupation :  0 
## insured_hobbies :  0 
## insured_relationship :  0 
## capital.gains :  0 
## capital.loss :  0 
## incident_date :  0 
## incident_type :  0 
## collision_type :  0 
## incident_severity :  0 
## authorities_contacted :  0 
## incident_state :  0 
## incident_city :  0 
## incident_location :  0 
## incident_hour_of_the_day :  0 
## number_of_vehicles_involved :  0 
## property_damage :  0 
## bodily_injuries :  0 
## witnesses :  0 
## police_report_available :  0 
## total_claim_amount :  0 
## injury_claim :  0 
## property_claim :  0 
## vehicle_claim :  0 
## auto_make :  0 
## auto_model :  0 
## auto_year :  0 
## fraud_reported :  0

c. Medidas Descriptivas

summary(df_aic)

##  months_as_customer      age        policy_number     policy_bind_date
##  Min.   :  0.0      Min.   :19.00   Min.   :100804   1/1/2006 :  3    
##  1st Qu.:115.8      1st Qu.:32.00   1st Qu.:335980   4/28/1992:  3    
##  Median :199.5      Median :38.00   Median :533135   8/5/1992 :  3    
##  Mean   :204.0      Mean   :38.95   Mean   :546239   1/28/2010:  2    
##  3rd Qu.:276.2      3rd Qu.:44.00   3rd Qu.:759100   1/29/1998:  2    
##  Max.   :479.0      Max.   :64.00   Max.   :999435   1/3/2004 :  2    
##                                                      (Other)  :985    
##  policy_state    policy_csl  policy_deductable policy_annual_premium
##  IL:338       100/300 :349   Min.   : 500      Min.   : 433.3       
##  IN:310       250/500 :351   1st Qu.: 500      1st Qu.:1089.6       
##  OH:352       500/1000:300   Median :1000      Median :1257.2       
##                              Mean   :1136      Mean   :1256.4       
##                              3rd Qu.:2000      3rd Qu.:1415.7       
##                              Max.   :2000      Max.   :2047.6       
##                                                                     
##  umbrella_limit      insured_zip     insured_sex  insured_education_level
##  Min.   :-1000000   Min.   :430104   FEMALE:537   Associate  :145        
##  1st Qu.:       0   1st Qu.:448404   MALE  :463   College    :122        
##  Median :       0   Median :466446                High School:160        
##  Mean   : 1101000   Mean   :501214                JD         :161        
##  3rd Qu.:       0   3rd Qu.:603251                Masters    :143        
##  Max.   :10000000   Max.   :620962                MD         :144        
##                                                   PhD        :125        
##          insured_occupation       insured_hobbies     insured_relationship
##  machine-op-inspct: 93      reading       : 64    husband       :170      
##  prof-specialty   : 85      exercise      : 57    not-in-family :174      
##  tech-support     : 78      paintball     : 57    other-relative:177      
##  exec-managerial  : 76      bungie-jumping: 56    own-child     :183      
##  sales            : 76      camping       : 55    unmarried     :141      
##  craft-repair     : 74      golf          : 55    wife          :155      
##  (Other)          :518      (Other)       :656                            
##  capital.gains     capital.loss       incident_date
##  Min.   :     0   Min.   :-111100   2/2/2015 : 28  
##  1st Qu.:     0   1st Qu.: -51500   2/17/2015: 26  
##  Median :     0   Median : -23250   1/7/2015 : 25  
##  Mean   : 25126   Mean   : -26794   1/10/2015: 24  
##  3rd Qu.: 51025   3rd Qu.:      0   1/24/2015: 24  
##  Max.   :100500   Max.   :      0   2/4/2015 : 24  
##                                     (Other)  :849  
##                   incident_type         collision_type      incident_severity
##  Multi-vehicle Collision :419   Front Collision:254    Major Damage  :276    
##  Parked Car              : 84   Rear Collision :470    Minor Damage  :354    
##  Single Vehicle Collision:403   Side Collision :276    Total Loss    :280    
##  Vehicle Theft           : 94                          Trivial Damage: 90    
##                                                                              
##                                                                              
##                                                                              
##  authorities_contacted incident_state     incident_city
##  Ambulance:196         NC:110         Arlington  :152  
##  Fire     :223         NY:262         Columbus   :149  
##  None     : 91         OH: 23         Hillsdale  :141  
##  Other    :198         PA: 30         Northbend  :145  
##  Police   :292         SC:248         Northbrook :122  
##                        VA:110         Riverwood  :134  
##                        WV:217         Springfield:157  
##         incident_location incident_hour_of_the_day number_of_vehicles_involved
##  1012 5th Lane   :  1     Min.   : 0.00            Min.   :1.000              
##  1028 Sky Lane   :  1     1st Qu.: 6.00            1st Qu.:1.000              
##  1030 Pine Lane  :  1     Median :12.00            Median :1.000              
##  1087 Flute Drive:  1     Mean   :11.64            Mean   :1.839              
##  1091 1st Drive  :  1     3rd Qu.:17.00            3rd Qu.:3.000              
##  1102 Apache Hwy :  1     Max.   :23.00            Max.   :4.000              
##  (Other)         :994                                                         
##  property_damage bodily_injuries   witnesses     police_report_available
##  NO :614         Min.   :0.000   Min.   :0.000   NO :598                
##  YES:386         1st Qu.:0.000   1st Qu.:1.000   YES:402                
##                  Median :1.000   Median :1.000                          
##                  Mean   :0.992   Mean   :1.487                          
##                  3rd Qu.:2.000   3rd Qu.:2.000                          
##                  Max.   :2.000   Max.   :3.000                          
##                                                                         
##  total_claim_amount  injury_claim   property_claim  vehicle_claim  
##  Min.   :   100     Min.   :    0   Min.   :    0   Min.   :   70  
##  1st Qu.: 41812     1st Qu.: 4295   1st Qu.: 4445   1st Qu.:30292  
##  Median : 58055     Median : 6775   Median : 6750   Median :42100  
##  Mean   : 52762     Mean   : 7433   Mean   : 7400   Mean   :37929  
##  3rd Qu.: 70592     3rd Qu.:11305   3rd Qu.:10885   3rd Qu.:50822  
##  Max.   :114920     Max.   :21450   Max.   :23670   Max.   :79560  
##                                                                    
##      auto_make      auto_model    auto_year    fraud_reported
##  Dodge    : 80   RAM     : 43   Min.   :1995   N:753         
##  Saab     : 80   Wrangler: 42   1st Qu.:2000   Y:247         
##  Suburu   : 80   A3      : 37   Median :2005                 
##  Nissan   : 78   Neon    : 37   Mean   :2005                 
##  Chevrolet: 76   MDX     : 36   3rd Qu.:2010                 
##  BMW      : 72   Jetta   : 35   Max.   :2015                 
##  (Other)  :534   (Other) :770

d. Medidas de Dispersión

medidas_dispersion <- data.frame(Variable = character(), Desviacion = double(), Varianza = double(), Rango = double(), IQR = double(), stringsAsFactors = FALSE)

for (nombre in names(df_aic)) {
  if (is.numeric(df_aic[[nombre]])) {
    # Calculo de medidas de dispersión
    desviacion <- formatC(sd(df_aic[[nombre]], na.rm = TRUE), format = "f", digits = 1)
    varianza <- formatC(var(df_aic[[nombre]], na.rm = TRUE), format = "f", digits = 1)
    rango <- formatC(max(df_aic[[nombre]], na.rm = TRUE) - min(df_aic[[nombre]], na.rm = TRUE), format = "f", digits = 1)
    iqr <- formatC(IQR(df_aic[[nombre]], na.rm = TRUE), format = "f", digits = 1)
    
    # Notación cientifica a max 2 decimales
    ajustarNotacion <- function(x) {
      numero <- as.numeric(x)
      if (abs(numero) > 1e+10) {
        return(formatC(numero, format = "e", digits = 1))
      } else {
        return(x)
      }
    }
    
    desviacion <- ajustarNotacion(desviacion)
    varianza <- ajustarNotacion(varianza)
    rango <- ajustarNotacion(rango)
    iqr <- ajustarNotacion(iqr)
    
    # Resultados a df
    medidas_dispersion <- rbind(medidas_dispersion, data.frame(Variable = nombre, Desviacion = desviacion, Varianza = varianza, Rango = rango, IQR = iqr, stringsAsFactors = FALSE))
  }
}

print(medidas_dispersion)

##                       Variable Desviacion     Varianza      Rango      IQR
## 1           months_as_customer      115.1      13251.0      479.0    160.5
## 2                          age        9.1         83.5       45.0     12.0
## 3                policy_number   257063.0      6.6e+10   898631.0 423119.5
## 4            policy_deductable      611.9     374378.4     1500.0   1500.0
## 5        policy_annual_premium      244.2      59617.7     1614.3    326.1
## 6               umbrella_limit  2297406.6      5.3e+12 11000000.0      0.0
## 7                  insured_zip    71701.6 5141121011.6   190858.0 154846.5
## 8                capital.gains    27872.2  776858847.6   100500.0  51025.0
## 9                 capital.loss    28104.1  789840250.6   111100.0  51500.0
## 10    incident_hour_of_the_day        7.0         48.3       23.0     11.0
## 11 number_of_vehicles_involved        1.0          1.0        3.0      2.0
## 12             bodily_injuries        0.8          0.7        2.0      2.0
## 13                   witnesses        1.1          1.2        3.0      1.0
## 14          total_claim_amount    26401.5  697040954.8   114820.0  28780.0
## 15                injury_claim     4881.0   23823691.0    21450.0   7010.0
## 16              property_claim     4824.7   23277982.7    23670.0   6440.0
## 17               vehicle_claim    18886.3  356690548.3    79490.0  20530.0
## 18                   auto_year        6.0         36.2       20.0     10.0

e. Transformación y Adicción de Variables

Columna de tipo de auto

# Cargar datos de tipos de vehículos
car_types <- read.csv("car_types.csv")

# Hacer merge con df_aic para añadir la columna vehicle_type
df_aic <- merge(df_aic, car_types, by = "auto_model", all.x = TRUE)

Transformar variables categóricas binarias

# Transformar 'police_report_available' de YES/NO a 1/0
df_aic$police_report_available <- ifelse(df_aic$police_report_available == "YES", 1, 0)

# Transformar 'property_damage' de YES/NO a 1/0
df_aic$property_damage <- ifelse(df_aic$property_damage == "YES", 1, 0)

# Transformar 'fraud_reported' de Y/N a 1/0
df_aic$fraud_reported <- ifelse(df_aic$fraud_reported == "Y", 1, 0)

# Transformar 'insured_sex' de MALE/FEMALE a 1/0
df_aic$insured_sex <- ifelse(df_aic$insured_sex == "MALE", 1, 0)

Dividir data de tiempo

df_aic <- df_aic %>%
  separate(incident_date, into = c("incident_day", "incident_month", "incident_year"), sep = "/")

Datos Geográficos de asegurado

# Base de datos geografícos por zip
uszips_data <- read_csv("uszips.csv")

# Preparación inicial
df_aic$insured_zip <- as.character(df_aic$insured_zip)
uszips_data$zip <- as.character(uszips_data$zip)
df_aic$modified_zip <- substr(df_aic$insured_zip, 1, nchar(df_aic$insured_zip) - 1)

# Añadir columnas vacías para los datos de uszips
df_aic$insured_lat <- NA
df_aic$insured_lng <- NA
df_aic$insured_county_name <- NA
df_aic$insured_city <- NA

# Función para intentar encontrar coincidencias ajustando el código postal
try_adjust_zip_and_merge <- function(row, uszips_data) {
  base_zip <- row$modified_zip
  for (i in c(-60:60)) {
    if (i == 0) next # Omitir el cero
    adjusted_zip <- as.character(as.numeric(base_zip) + i)
    
    # Buscar coincidencia en uszips_data
    match <- uszips_data[uszips_data$zip == adjusted_zip, ]
    if (nrow(match) > 0) {
      row$insured_lat <- match$lat[1]
      row$insured_lng <- match$lng[1]
      row$insured_county_name <- match$county_name[1]
      row$insured_city <- match$city[1]
      break # Salir del bucle si se encuentra una coincidencia
    }
  }
  return(row)
}

# Aplicar la función a cada fila de df_aic
for (i in 1:nrow(df_aic)) {
  df_aic[i, ] <- try_adjust_zip_and_merge(df_aic[i, ], uszips_data)
}

# Limpiar el dataframe final
df_aic$modified_zip <- NULL

f. Creación de Dataframe Espacial

# Opciones para tigris
options(tigris_use_cache = TRUE, tigris_class = "sf")

# Obtener datos geográficos para condados en Ohio, Illinois, e Indiana
states_of_interest <- c("Ohio", "Illinois", "Indiana")
areas_sf <- counties(state = states_of_interest, cb = TRUE)

## Retrieving data for the year 2022

# Función para calcular la moda
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

# Crear un nuevo dataframe con promedios agrupados por insured_county_name
df_aic_med <- df_aic %>%
  group_by(insured_county_name) %>%
  summarise(
    med_total_claim_amount = median(total_claim_amount, na.rm = TRUE),
    med_policy_annual_premium = median(policy_annual_premium, na.rm = TRUE),
    med_umbrella_limit = median(umbrella_limit, na.rm = TRUE),
    med_vehicle_claim = median(vehicle_claim, na.rm = TRUE),
    med_injury_claim = median(injury_claim + 0.01, na.rm = TRUE),
    med_number_of_vehicles_involved = median(number_of_vehicles_involved, na.rm = TRUE),
    med_incident_hour_of_the_day = median(incident_hour_of_the_day + 0.01, na.rm = TRUE), 
    # Variables categoricas
    mode_police_report_available = Mode(police_report_available),
    mode_property_damage = Mode(property_damage)
  )

# Realizar el join entre areas_sf y df_aic_med
df_aic_sf <- areas_sf %>%
  left_join(df_aic_med, by = c("NAME" = "insured_county_name"))

# Configurar tmap para visualización estática
tmap_mode("plot")

## tmap mode set to plotting

g. Identificación de Patrones

y/o tendencias en los datos mediante el uso de gráficos incluyendo bar plots, line plots, pie plots, histogramas, matriz de correlación, box plot, scatter plot, qq- plot, etc Mostrar al menos 4 – 6 gráficos.

Distribución de Variables Cuantitativas

# Identificar las columnas numéricas en df_aic
columnas_numericas <- sapply(df_aic, is.numeric)

# Crear histogramas para cada columna numérica
for (nombre_columna in names(columnas_numericas[columnas_numericas])) {
  p <- ggplot(df_aic, aes_string(x = nombre_columna)) +
    geom_histogram(bins = 30, fill = "blue", color = "black") + # Ajusta los bins según necesites
    theme_minimal() +
    labs(title = paste("Histogram of", nombre_columna),
         x = nombre_columna,
         y = "Frequency")
  
  print(p) # Mostrar el histograma
}

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Boxplot de Reclamo Total vs Tipo de Carro

ggplot(df_aic, aes(x = type, y = total_claim_amount, fill = type)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Diagrama de caja del Monto Total de la Reclamación por Tipo",
       x = "Tipo de Auto",
       y = "Monto Total de Reclamo") + 
  scale_fill_brewer(palette = "Pastel1")

Histograma del Monto Total del Reclamo con Indicador de Daños a la Propiedad

ggplot(df_aic, aes(x = total_claim_amount, fill = as.factor(property_damage))) +
  geom_histogram(position = "dodge", binwidth = 5000) + 
  scale_fill_manual(values = c("0" = "navy", "1" = "lightblue")) +
  labs(title = "Histograma del Monto Total del Reclamo con Indicador de Daños a la Propiedad",
       x = "Monto Total de Reclamo",
       y = "Numero de casos",
       fill = "Daños a la Propiedad") +
  theme_minimal()

Gráfico de Dispersión por Reclamo Total vs Reclamo por Accidente

# Crear scatter plot con ggplot2 utilizando una escala de color de azul claro a azul marino por age
ggplot(df_aic, aes(x = total_claim_amount, y = injury_claim, color = age)) +
  geom_point() + # Añadir puntos al scatter plot
  theme_minimal() + # Aplicar un tema minimalista
  scale_color_gradient(low = "lightblue", high = "navy") + # Gradiente de color de azul claro a azul marino
  labs(title = "Monto Total del Reclamo frente al Reclamo por Accidente por Edad",
       x = "Monto total del Reclamo",
       y = "Reclamación por Accidente",
       color = "Edad") # Personalizar las etiquetas

Box Plot for Number of Vehicles Involved

ggplot(df_aic, aes(x = as.factor(number_of_vehicles_involved), y = total_claim_amount, fill = as.factor(number_of_vehicles_involved))) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Monto Total del Reclamo por Número de Vehículos Involucrados",
       x = "Numero de Vehículos Involucrados",
       y = "Monto Total de Reclamo") +
  scale_x_discrete(name = "Numero de Vehículos Involucrados") +
  scale_fill_brewer(palette = "Pastel1", name = "Numero de Vehículos Involucrados") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")  

Matriz de correlación

# Crear df_corrplot con solo variables numéricas de df_aic
df_corrplot <- df_aic[sapply(df_aic, is.numeric)]

# Calcular la matriz de correlación
matriz_corr <- cor(df_corrplot, use = "complete.obs") # 'use' maneja los valores NA

# Generar el gráfico de correlación
corrplot(matriz_corr, order = "AOE", tl.cex = 0.7, tl.col = "black") # Ajusta los parámetros según necesidad

Monto de Accidente y Poliza Annual por Condado

# Mapa para Avg Total Claim Amount
total_claim_map <- tm_shape(df_aic_sf) + 
  tm_fill("med_total_claim_amount", palette = "Blues", style = "quantile", title = "Monto Total Promedio del Reclamo") +
  tm_borders(alpha = .4) + 
  tm_layout(legend.text.size = 0.8, legend.title.size = 1.1, frame = FALSE)

# Mapa para Avg Policy Annual Premium
annual_premium_map <- tm_shape(df_aic_sf) + 
  tm_fill("med_policy_annual_premium", palette = "BuPu", style = "quantile", title = "Prima Anual Promedio de la Póliza") +
  tm_borders(alpha = .4) + 
  tm_layout(legend.text.size = 0.8, legend.title.size = 1.1, frame = FALSE)

# Arreglar mapas lado a lado
tmap_arrange(total_claim_map, annual_premium_map, nrow = 1)

## Some legend labels were too wide. These labels have been resized to 0.73. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Umbrella Limit per State

# Ajusta estos valores según cómo estén representados los estados en tus datos
states_codes <- c("Ohio" = "39", "Illinois" = "17", "Indiana" = "18")

# Crear un mapa para cada estado
for(state_name in names(states_codes)) {
  state_code <- states_codes[state_name]
  
  # Filtrar df_aic_sf por estado actual en el bucle
  state_data <- df_aic_sf %>% 
    filter(STATEFP == state_code)  
  
  # Crear el mapa con la leyenda en la esquina superior derecha de la imagen
  map <- tm_shape(state_data) + 
    tm_polygons(col = "med_umbrella_limit", 
                title = paste("Umbrella Limit in", state_name)) + 
    tm_borders() +
    tm_layout(
      frame = FALSE,
      title = paste(""),
      legend.position = c("right", "top"),
      legend.just = c(1, 1),  # Ajusta la leyenda para que el punto de referencia sea la esquina superior derecha
      legend.outside = TRUE,  # Coloca la leyenda fuera del área del mapa
      legend.outside.position = "right",  # Coloca la leyenda en el lado derecho fuera del área del mapa
      outer.margins = c(0, 0, 0, 0),  # Ajusta si es necesario para los márgenes exteriores
      legend.frame = TRUE
    )
  
  # Imprimir el mapa
  print(map)
}

Especificación Preferida

A partir de los resultados de EDA describir la especificación del modelo de regresión lineal a estimar. Brevemente, describir cómo es el posible impacto de cada una de las variables explicativas sobre la principal variable de estudio.

Estimación y Modelado

Estimación de cada uno de los siguientes modelos de Supervised Machine Learning (SML): ## a. OLS Regresión ### Modelo

model_ols <- lm(total_claim_amount ~ police_report_available + property_damage + log(vehicle_claim) + log(injury_claim + 0.01) + number_of_vehicles_involved + incident_hour_of_the_day, data = df_aic)

# Mostrar el resumen del modelo
summary(model_ols)

## 
## Call:
## lm(formula = total_claim_amount ~ police_report_available + property_damage + 
##     log(vehicle_claim) + log(injury_claim + 0.01) + number_of_vehicles_involved + 
##     incident_hour_of_the_day, data = df_aic)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -22312  -6925  -1140   5905  94162 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 -197430.70    3310.75 -59.633  < 2e-16 ***
## police_report_available        -419.40     651.75  -0.643 0.520054    
## property_damage                1086.78     656.17   1.656 0.097989 .  
## log(vehicle_claim)            23864.53     374.91  63.654  < 2e-16 ***
## log(injury_claim + 0.01)        975.69     149.14   6.542 9.71e-11 ***
## number_of_vehicles_involved   -1152.10     329.70  -3.494 0.000496 ***
## incident_hour_of_the_day        -40.31      46.97  -0.858 0.390989    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9950 on 993 degrees of freedom
## Multiple R-squared:  0.8588, Adjusted R-squared:  0.858 
## F-statistic:  1007 on 6 and 993 DF,  p-value: < 2.2e-16

RMSE

df_aic_test <- df_aic

# Asegúrate de que las transformaciones se aplican correctamente
df_aic_test$log_vehicle_claim <- log(abs(df_aic$vehicle_claim))
df_aic_test$log_injury_claim_plus <- log(abs(df_aic$injury_claim + 0.01))

# Ahora, ajusta el modelo con las variables transformadas correctamente
model_ols <- lm(total_claim_amount ~ police_report_available + property_damage + log_vehicle_claim + log_injury_claim_plus + number_of_vehicles_involved + incident_hour_of_the_day, data = df_aic_test)

# Calcula las predicciones
predicciones <- predict(model_ols, df_aic_test)

# Calcula los residuos y el RMSE como antes
residuos <- df_aic_test$total_claim_amount - predicciones
rmse <- sqrt(mean(residuos^2))

print(rmse)

## [1] 9914.873

b. SAR

Modelo

df_aic_sf <- na.omit(df_aic_sf)
df_aic_sp <- as(df_aic_sf, "Spatial")
nb <- poly2nb(df_aic_sp)
listw <- nb2listw(nb, style="W", zero.policy = TRUE)

model_sar <- lagsarlm(med_total_claim_amount ~ mode_police_report_available + mode_property_damage +
                      log(med_vehicle_claim + 0.01) + log(med_injury_claim + 0.01) +
                      med_number_of_vehicles_involved + med_incident_hour_of_the_day,
                      data = df_aic_sf, listw = listw, zero.policy = TRUE, method="eigen")

# Mostrar el resumen del modelo
summary(model_sar)

## 
## Call:
## lagsarlm(formula = med_total_claim_amount ~ mode_police_report_available + 
##     mode_property_damage + log(med_vehicle_claim + 0.01) + log(med_injury_claim + 
##     0.01) + med_number_of_vehicles_involved + med_incident_hour_of_the_day, 
##     data = df_aic_sf, listw = listw, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -12954.85  -4321.42   -846.09   4095.68  29493.97 
## 
## Type: lag 
## Regions with no neighbours included:
##  76 
## Coefficients: (numerical Hessian approximate standard errors) 
##                                    Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept)                     -229677.541    8815.861 -26.0528 < 2.2e-16
## mode_police_report_available        162.126     929.008   0.1745    0.8615
## mode_property_damage               4168.504     957.471   4.3537 1.339e-05
## log(med_vehicle_claim + 0.01)     26896.088     982.178  27.3841 < 2.2e-16
## log(med_injury_claim + 0.01)        327.537     498.252   0.6574    0.5109
## med_number_of_vehicles_involved     536.047     531.847   1.0079    0.3135
## med_incident_hour_of_the_day       -136.929      94.659  -1.4465    0.1480
## 
## Rho: -0.020858, LR test value: 0.19101, p-value: 0.66207
## Approximate (numerical Hessian) standard error: 0.049067
##     z-value: -0.42508, p-value: 0.67078
## Wald statistic: 0.1807, p-value: 0.67078
## 
## Log likelihood: -2440.081 for lag model
## ML residual variance (sigma squared): 43190000, (sigma: 6571.9)
## Number of observations: 239 
## Number of parameters estimated: 9 
## AIC: 4898.2, (AIC for lm: 4896.4)

c. SEM

Modelo

model_sem <- errorsarlm(med_total_claim_amount ~ mode_police_report_available + mode_property_damage +
                        log(med_vehicle_claim + 0.01) + log(med_injury_claim + 0.01) +
                        med_number_of_vehicles_involved + med_incident_hour_of_the_day,
                        data = df_aic_sf, listw = listw, zero.policy = TRUE, method="eigen")

# Mostrar el resumen del modelo
summary(model_sem)

## 
## Call:
## errorsarlm(formula = med_total_claim_amount ~ mode_police_report_available + 
##     mode_property_damage + log(med_vehicle_claim + 0.01) + log(med_injury_claim + 
##     0.01) + med_number_of_vehicles_involved + med_incident_hour_of_the_day, 
##     data = df_aic_sf, listw = listw, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -12746.94  -4361.60   -965.61   4176.57  29301.70 
## 
## Type: error 
## Regions with no neighbours included:
##  76 
## Coefficients: (asymptotic standard errors) 
##                                    Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept)                     -229367.767    8754.928 -26.1987 < 2.2e-16
## mode_police_report_available        140.456     922.599   0.1522    0.8790
## mode_property_damage               4002.041     949.345   4.2156 2.491e-05
## log(med_vehicle_claim + 0.01)     26724.135     967.805  27.6131 < 2.2e-16
## log(med_injury_claim + 0.01)        330.366     486.831   0.6786    0.4974
## med_number_of_vehicles_involved     644.446     518.224   1.2436    0.2137
## med_incident_hour_of_the_day       -117.434      93.169  -1.2604    0.2075
## 
## Lambda: 0.12704, LR test value: 1.6709, p-value: 0.19614
## Asymptotic standard error: 0.092272
##     z-value: 1.3768, p-value: 0.16859
## Wald statistic: 1.8955, p-value: 0.16859
## 
## Log likelihood: -2439.341 for error model
## ML residual variance (sigma squared): 42775000, (sigma: 6540.2)
## Number of observations: 239 
## Number of parameters estimated: 9 
## AIC: 4896.7, (AIC for lm: 4896.4)

d. XGBoost Regresión

Modelo

# Transformaciones Necesarias
df_aic_transformed <- df_aic %>%
  mutate(
    log_vehicle_claim = log(vehicle_claim),
    log_injury_claim = log(injury_claim + 0.01) 
  ) 

df_aic_transformed <- df_aic_transformed[c("total_claim_amount", "police_report_available", "property_damage", "log_vehicle_claim", "log_injury_claim", "number_of_vehicles_involved", "incident_hour_of_the_day")]


set.seed(123)
cv_data <- createDataPartition(y = df_aic_transformed$total_claim_amount, p = 0.7, list = FALSE)
cv_train <- df_aic_transformed[cv_data, ]
cv_test <- df_aic_transformed[-cv_data, ]

# Preparación de Matrices
train_x <- data.matrix(cv_train[, -1])
train_y <- cv_train[, 1]
test_x <- data.matrix(cv_test[, -1])
test_y <- cv_test[, 1]

# Dataframes de prueba y test
xgb_train <- xgb.DMatrix(data = train_x, label = train_y)
xgb_test <- xgb.DMatrix(data = test_x, label = test_y)

# XGBoost regression model
watchlist <- list(train = xgb_train, test = xgb_test)
model_xgb <- xgb.train(data = xgb_train, objective = "reg:squarederror", max.depth = 3, watchlist = watchlist, nrounds = 70, eta = 0.1)

## [1]  train-rmse:53326.710755 test-rmse:52913.885784 
## [2]  train-rmse:48107.326928 test-rmse:47728.422530 
## [3]  train-rmse:43411.039591 test-rmse:43086.550037 
## [4]  train-rmse:39179.303365 test-rmse:38898.496021 
## [5]  train-rmse:35372.078861 test-rmse:35143.585459 
## [6]  train-rmse:31944.359781 test-rmse:31734.675487 
## [7]  train-rmse:28860.347040 test-rmse:28710.842077 
## [8]  train-rmse:26082.581516 test-rmse:25945.354058 
## [9]  train-rmse:23582.811195 test-rmse:23511.432088 
## [10] train-rmse:21335.623534 test-rmse:21318.688142 
## [11] train-rmse:19313.266654 test-rmse:19318.124343 
## [12] train-rmse:17494.924734 test-rmse:17462.527193 
## [13] train-rmse:15859.764347 test-rmse:15853.962855 
## [14] train-rmse:14387.767287 test-rmse:14425.427268 
## [15] train-rmse:13068.264632 test-rmse:13144.498498 
## [16] train-rmse:11887.140539 test-rmse:11959.607085 
## [17] train-rmse:10828.483662 test-rmse:10896.316712 
## [18] train-rmse:9872.921658  test-rmse:9965.033993 
## [19] train-rmse:9024.341545  test-rmse:9152.524432 
## [20] train-rmse:8263.612845  test-rmse:8396.536582 
## [21] train-rmse:7590.252130  test-rmse:7748.553680 
## [22] train-rmse:6993.412375  test-rmse:7161.530409 
## [23] train-rmse:6454.445322  test-rmse:6639.241159 
## [24] train-rmse:5977.644473  test-rmse:6176.509372 
## [25] train-rmse:5556.707861  test-rmse:5773.714910 
## [26] train-rmse:5190.004040  test-rmse:5411.597963 
## [27] train-rmse:4865.921319  test-rmse:5110.647581 
## [28] train-rmse:4588.944990  test-rmse:4852.145076 
## [29] train-rmse:4343.886103  test-rmse:4621.515823 
## [30] train-rmse:4134.441112  test-rmse:4413.928111 
## [31] train-rmse:3951.595868  test-rmse:4249.661727 
## [32] train-rmse:3797.695716  test-rmse:4106.489675 
## [33] train-rmse:3663.356204  test-rmse:3985.738501 
## [34] train-rmse:3551.297994  test-rmse:3878.377582 
## [35] train-rmse:3453.834316  test-rmse:3790.878192 
## [36] train-rmse:3369.821607  test-rmse:3715.278973 
## [37] train-rmse:3297.281593  test-rmse:3642.737500 
## [38] train-rmse:3237.014889  test-rmse:3587.086591 
## [39] train-rmse:3181.347440  test-rmse:3533.699532 
## [40] train-rmse:3136.168693  test-rmse:3495.784510 
## [41] train-rmse:3094.416063  test-rmse:3458.343779 
## [42] train-rmse:3062.625207  test-rmse:3429.142526 
## [43] train-rmse:3031.545871  test-rmse:3401.522322 
## [44] train-rmse:3003.830008  test-rmse:3379.350245 
## [45] train-rmse:2984.901069  test-rmse:3364.762494 
## [46] train-rmse:2963.032642  test-rmse:3348.920911 
## [47] train-rmse:2944.198309  test-rmse:3334.213945 
## [48] train-rmse:2930.070104  test-rmse:3323.911990 
## [49] train-rmse:2913.435116  test-rmse:3317.362620 
## [50] train-rmse:2903.450271  test-rmse:3308.855638 
## [51] train-rmse:2894.157913  test-rmse:3304.039850 
## [52] train-rmse:2880.487650  test-rmse:3300.440099 
## [53] train-rmse:2869.475342  test-rmse:3295.016747 
## [54] train-rmse:2862.072531  test-rmse:3289.625826 
## [55] train-rmse:2854.864375  test-rmse:3286.107156 
## [56] train-rmse:2846.987900  test-rmse:3283.997066 
## [57] train-rmse:2842.064021  test-rmse:3282.171625 
## [58] train-rmse:2832.984051  test-rmse:3280.796569 
## [59] train-rmse:2828.535467  test-rmse:3279.886700 
## [60] train-rmse:2822.185214  test-rmse:3277.683291 
## [61] train-rmse:2818.744740  test-rmse:3276.191652 
## [62] train-rmse:2815.384078  test-rmse:3275.836318 
## [63] train-rmse:2805.038377  test-rmse:3272.029084 
## [64] train-rmse:2799.497590  test-rmse:3271.795408 
## [65] train-rmse:2793.761314  test-rmse:3271.930104 
## [66] train-rmse:2783.824131  test-rmse:3261.833695 
## [67] train-rmse:2776.416844  test-rmse:3265.465576 
## [68] train-rmse:2773.288009  test-rmse:3266.782008 
## [69] train-rmse:2769.526196  test-rmse:3267.147256 
## [70] train-rmse:2765.288923  test-rmse:3266.365310

# Estimar modelo final
reg_xgb <- xgboost(data = xgb_train, objective = "reg:squarederror", max.depth = 3, nrounds = 59, verbose = 0)

# Make predictions and calculate RMSE
prediction_xgb_test <- predict(reg_xgb, xgb_test)
RMSE_XGB <- sqrt(mean((prediction_xgb_test - test_y)^2))
print(paste("RMSE XGBoost:", RMSE_XGB))

## [1] "RMSE XGBoost: 3339.40193750338"

# Ensure predictions and test labels are finite
if(all(is.finite(prediction_xgb_test)) && all(is.finite(test_y))) {
  # Diagnostic checks of regression residuals
  xgb_reg_residuals <- test_y - prediction_xgb_test
  
  # Check if residuals contain only finite values before plotting
  if(all(is.finite(xgb_reg_residuals))) {
    plot(xgb_reg_residuals, xlab = "Predicted Values", ylab = "Residuals", main = "XGBoost Regression Residuals")
    abline(h = 0, col = "red") # Ensure horizontal line at 0
  } else {
    cat("Residuals contain non-finite values; cannot plot.\n")
  }
} else {
  cat("Predictions or test labels contain non-finite values; cannot calculate RMSE or plot residuals.\n")
}

# Plot the importance of variables
importance_matrix <- xgb.importance(model = reg_xgb)
xgb.plot.importance(importance_matrix)

RMSE

print(RMSE_XGB)

## [1] 3339.402

e. Decision Trees

Modelo

model_dt <- rpart(total_claim_amount ~ police_report_available + property_damage + 
                  log(vehicle_claim) + log(injury_claim + 0.01) + 
                  number_of_vehicles_involved + incident_hour_of_the_day, 
                  data = df_aic, 
                  method = "anova") # Usar "anova" para regresión

# Mostrar el resumen del modelo
summary(model_dt)

## Call:
## rpart(formula = total_claim_amount ~ police_report_available + 
##     property_damage + log(vehicle_claim) + log(injury_claim + 
##     0.01) + number_of_vehicles_involved + incident_hour_of_the_day, 
##     data = df_aic, method = "anova")
##   n= 1000 
## 
##           CP nsplit  rel error     xerror        xstd
## 1 0.71980532      0 1.00000000 1.00128183 0.039313456
## 2 0.14041228      1 0.28019468 0.28449191 0.013292832
## 3 0.03203047      2 0.13978240 0.14962416 0.007622793
## 4 0.02814189      3 0.10775193 0.10984892 0.006229195
## 5 0.02532398      4 0.07961004 0.08619526 0.005439850
## 6 0.01000000      5 0.05428606 0.06052631 0.003335295
## 
## Variable importance
##          log(vehicle_claim)    log(injury_claim + 0.01) 
##                          58                          41 
## number_of_vehicles_involved 
##                           1 
## 
## Node number 1: 1000 observations,    complexity param=0.7198053
##   mean=52761.94, MSE=6.963439e+08 
##   left son=2 (216 obs) right son=3 (784 obs)
##   Primary splits:
##       log(vehicle_claim)          < 10.17656 to the left,  improve=0.71980530, (0 missing)
##       log(injury_claim + 0.01)    < 8.25449  to the left,  improve=0.66695230, (0 missing)
##       incident_hour_of_the_day    < 9.5      to the left,  improve=0.09858620, (0 missing)
##       number_of_vehicles_involved < 1.5      to the left,  improve=0.08157524, (0 missing)
##       property_damage             < 0.5      to the left,  improve=0.01552611, (0 missing)
##   Surrogate splits:
##       log(injury_claim + 0.01) < 8.161948 to the left,  agree=0.957, adj=0.801, (0 split)
## 
## Node number 2: 216 observations,    complexity param=0.03203047
##   mean=10108.84, MSE=1.102546e+08 
##   left son=4 (180 obs) right son=5 (36 obs)
##   Primary splits:
##       log(vehicle_claim)          < 9.626001 to the left,  improve=0.93656210, (0 missing)
##       log(injury_claim + 0.01)    < 7.620277 to the left,  improve=0.91224760, (0 missing)
##       number_of_vehicles_involved < 2        to the left,  improve=0.51995230, (0 missing)
##       incident_hour_of_the_day    < 9.5      to the left,  improve=0.14712030, (0 missing)
##       police_report_available     < 0.5      to the left,  improve=0.01810053, (0 missing)
##   Surrogate splits:
##       log(injury_claim + 0.01)    < 7.620277 to the left,  agree=0.986, adj=0.917, (0 split)
##       number_of_vehicles_involved < 2        to the left,  agree=0.931, adj=0.583, (0 split)
## 
## Node number 3: 784 observations,    complexity param=0.1404123
##   mean=64513.3, MSE=2.184909e+08 
##   left son=6 (469 obs) right son=7 (315 obs)
##   Primary splits:
##       log(vehicle_claim)          < 10.78457 to the left,  improve=0.570794100, (0 missing)
##       log(injury_claim + 0.01)    < 9.472321 to the left,  improve=0.326795200, (0 missing)
##       incident_hour_of_the_day    < 0.5      to the left,  improve=0.004354350, (0 missing)
##       number_of_vehicles_involved < 1.5      to the right, improve=0.004039198, (0 missing)
##       property_damage             < 0.5      to the left,  improve=0.003338463, (0 missing)
##   Surrogate splits:
##       log(injury_claim + 0.01) < 9.413282 to the left,  agree=0.724, adj=0.314, (0 split)
## 
## Node number 4: 180 observations
##   mean=5564.389, MSE=4814187 
## 
## Node number 5: 36 observations
##   mean=32831.11, MSE=1.789498e+07 
## 
## Node number 6: 469 observations,    complexity param=0.02814189
##   mean=55361.11, MSE=7.990679e+07 
##   left son=12 (185 obs) right son=13 (284 obs)
##   Primary splits:
##       log(vehicle_claim)          < 10.55229 to the left,  improve=0.522902200, (0 missing)
##       log(injury_claim + 0.01)    < 9.23308  to the left,  improve=0.293596700, (0 missing)
##       incident_hour_of_the_day    < 20.5     to the left,  improve=0.008676308, (0 missing)
##       property_damage             < 0.5      to the right, improve=0.006231092, (0 missing)
##       number_of_vehicles_involved < 3.5      to the left,  improve=0.001643506, (0 missing)
##   Surrogate splits:
##       log(injury_claim + 0.01) < 8.523175 to the left,  agree=0.684, adj=0.2, (0 split)
## 
## Node number 7: 315 observations,    complexity param=0.02532398
##   mean=78139.9, MSE=1.144297e+08 
##   left son=14 (201 obs) right son=15 (114 obs)
##   Primary splits:
##       log(vehicle_claim)          < 10.96136 to the left,  improve=0.4892227000, (0 missing)
##       log(injury_claim + 0.01)    < 9.574289 to the left,  improve=0.2729756000, (0 missing)
##       number_of_vehicles_involved < 1.5      to the right, improve=0.0162950400, (0 missing)
##       incident_hour_of_the_day    < 20.5     to the left,  improve=0.0059580980, (0 missing)
##       property_damage             < 0.5      to the left,  improve=0.0007988186, (0 missing)
##   Surrogate splits:
##       log(injury_claim + 0.01) < 9.686568 to the left,  agree=0.692, adj=0.149, (0 split)
## 
## Node number 12: 185 observations
##   mean=47352.16, MSE=3.459087e+07 
## 
## Node number 13: 284 observations
##   mean=60578.2, MSE=4.042445e+07 
## 
## Node number 14: 201 observations
##   mean=72505.12, MSE=4.294288e+07 
## 
## Node number 15: 114 observations
##   mean=88074.91, MSE=8.578616e+07

plot(model_dt, uniform = TRUE, main = "Árbol de Decisión")
text(model_dt, use.n = TRUE)

RMSE

# Ajuste del modelo, asegurándose de que las transformaciones están aplicadas correctamente
df_aic$log_vehicle_claim <- log(df_aic$vehicle_claim)
df_aic$log_injury_claim_plus <- log(df_aic$injury_claim + 0.01)

model_dt <- rpart(total_claim_amount ~ police_report_available + property_damage + 
                  log_vehicle_claim + log_injury_claim_plus + 
                  number_of_vehicles_involved + incident_hour_of_the_day, 
                  data = df_aic, 
                  method = "anova")

# Predicciones con el modelo de árbol de decisión, utilizando las mismas transformaciones
predicciones_dt <- predict(model_dt, df_aic)

# Calcular los residuos
residuos_dt <- df_aic$total_claim_amount - predicciones_dt

# Calcular el RMSE
rmse_dt <- sqrt(mean(residuos_dt^2))

# Imprimir el RMSE
print(rmse_dt)

## [1] 6148.314

f. Random Forest

Modelo

model_rf <- randomForest(total_claim_amount ~ police_report_available + property_damage + vehicle_claim + injury_claim + number_of_vehicles_involved + incident_hour_of_the_day, data= df_aic, proximity=TRUE)

# Mostrar el resumen del modelo
print(model_rf)

## 
## Call:
##  randomForest(formula = total_claim_amount ~ police_report_available +      property_damage + vehicle_claim + injury_claim + number_of_vehicles_involved +      incident_hour_of_the_day, data = df_aic, proximity = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 17259368
##                     % Var explained: 97.52

# Visualizar la importancia de las variables
varImpPlot(model_rf)

### RMSE

# Hacer predicciones con el modelo de Random Forest
predicciones_rf <- predict(model_rf, df_aic)

# Calcular los residuos
residuos_rf <- df_aic$total_claim_amount - predicciones_rf

# Calcular el RMSE
rmse_rf <- sqrt(mean(residuos_rf^2))

# Imprimir el RMSE
print(rmse_rf)

## [1] 2373.398

g. Neural Networks Regresión

df_aic$vehicle_claim <- scale(df_aic$vehicle_claim)
df_aic$injury_claim <- scale(df_aic$injury_claim)

nn_model <- neuralnet(total_claim_amount ~ police_report_available + property_damage + 
                      vehicle_claim + injury_claim + 
                      number_of_vehicles_involved + incident_hour_of_the_day, 
                      data = df_aic, 
                      hidden = c(5, 3), 
                      linear.output = TRUE,
                      stepmax = 1000000) # Valor incrementado

# Intenta visualizar nuevamente
plot(nn_model)

Pruebas de Diagnóstico

Pruebas de Diagnóstico de los Resultados Obtenidos de la Estmación de Modelos de Regresión

a. Multicolinealidad

vif_resultados <- vif(model_ols)

# Mostrar los resultados del VIF
print(vif_resultados)

##     police_report_available             property_damage 
##                    1.031500                    1.030786 
##           log_vehicle_claim       log_injury_claim_plus 
##                    1.419454                    1.198633 
## number_of_vehicles_involved    incident_hour_of_the_day 
##                    1.138766                    1.075865

b. Heterocedasticidad

bp_test_resultado <- bptest(model_ols)

# Mostrar los resultados del test
print(bp_test_resultado)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 16.109, df = 6, p-value = 0.01318

c. Autocorrelación Serial

NA

d. Autocorrelación Espacial

residuos_sar <- residuals(model_sar)

# Realizar el test de Moran para los residuos
moran_test_resultado <- moran.test(residuos_sar, listw)

# Mostrar los resultados del test de Moran
print(moran_test_resultado)

## 
##  Moran I test under randomisation
## 
## data:  residuos_sar  
## weights: listw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.614, p-value = 0.05327
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.064915215      -0.004219409       0.001834886

e. Normalidad de los Residuales

residuos <- residuals(model_ols)

# Realizar el test de Shapiro-Wilk para normalidad
shapiro_test <- shapiro.test(residuos)

# Mostrar los resultados del test
print(shapiro_test)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.93996, p-value < 2.2e-16

Nota: En caso de que las pruebas de diagnóstico identifiquen cualquiera de los anteriores a) – e) plantear una solución para mejorar la estimación de la especificiación del modelo.

Multicolinealidad

• Los resultados muestran que no existe multicolinealidad en el modelo, lo cual es una buena señal. Los valores de VIF para las variables police_report_available, property_damage, log(vehicle_claim), log(injury_claim + 0.01), number_of_vehicles_involved, e incident_hour_of_the_day están todos por debajo del umbral comúnmente aceptado de 5 o 10, indicando que no hay preocupaciones significativas de multicolinealidad entre estas variables.

Heteroscedasticidad

• Se identificó la presencia de heteroscedasticidad en el modelo a través del test de Breusch-Pagan, lo que sugiere que la varianza de los errores no es constante. Una posible solución para mitigar este problema es transformar las variables dependientes utilizando logaritmos, lo que puede ayudar a estabilizar la varianza a lo largo de los datos. Otra estrategia podría ser el uso de técnicas de ponderación o la aplicación de modelos robustos que sean menos sensibles a la heteroscedasticidad.

Ausencia de Necesidad de Contemplar Autocorrelación Serial

• Dado que la base de datos es de corte transversal y no utiliza series de tiempo, no es necesario considerar la autocorrelación serial en el análisis. Esto simplifica el análisis ya que la autocorrelación serial puede introducir sesgos significativos en los modelos de series de tiempo.

Autocorrelación Espacial

• Se exploró la autocorrelación espacial mediante modelos SEM y SAR, utilizando el test de Moran para evaluar el impacto espacial en los residuos. Los resultados indican que, a primera vista, no hay un impacto espacial significativo (p-valor = 0.05327), lo cual sugiere que la estructura espacial de los datos no influye de manera considerable en el modelo.

Normalidad de los Residuos

• El test de Shapiro-Wilk revela que los residuos del modelo no se distribuyen normalmente (p-valor = 2.544e-10). Esto podría ser problemático para ciertas inferencias estadísticas que asumen normalidad de los residuos. Como posibles soluciones, se podrían considerar transformaciones de los datos para mejorar la normalidad, el uso de métodos no paramétricos, o ajustar el modelo para manejar mejor las distribuciones de los errores no normales.

Evaluación y Selección de Modelo de Regresión

Comparación de RMSE

Presentar los valores de la métrica RMSE de cada uno de los modelos estimados en 4) en un gráfico de barras.

Modelo	RMSE
Random Forest	47558.14
Decision Tree	18211.22
XGBoost	3339.402
OLS	23931.76

datos <- data.frame(
  Modelo = c("Random Forest", "Decision Tree", "XGBoost", "OLS"),
  RMSE = c(47558.14, 18211.22, 3339.402, 23931.76)
)

# Generar el gráfico de barras
ggplot(datos, aes(x = Modelo, y = RMSE, fill = Modelo)) +
  geom_bar(stat = "identity", show.legend = FALSE) +
  theme_minimal() +
  labs(title = "RMSE por Modelo",
       x = "Modelo",
       y = "RMSE") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Me decantaría por el modelo XGBoost por las siguientes razones concretas:

• Menor RMSE: XGBoost lidera con el menor RMSE, señal de que sus predicciones suelen estar más cerca de los valores reales, un factor crítico para la confiabilidad en un modelo de predicción de reclamaciones.

• Habilidad para capturar complejidades: Se destaca en descifrar interacciones entre variables que otros modelos podrían no captar tan eficientemente, lo que es clave en un contexto donde las relaciones causales pueden ser intrincadas.

• Regularización incorporada: La regularización es una salvaguarda contra el sobreajuste. Esto significa que XGBoost no solo se ajusta bien a los datos con los que se entrena, sino que también mantiene su rendimiento en datos no vistos, lo cual es vital para aplicaciones prácticas.

• Flexibilidad en la configuración: Ofrece una variedad de parámetros que pueden ajustarse para mejorar el rendimiento y adaptarse mejor a las características específicas del conjunto de datos.

Dicho esto, es importante recordar que XGBoost puede sacrificar interpretabilidad por rendimiento; los modelos basados en árboles son a menudo menos transparentes que, por ejemplo, un modelo lineal. Por lo tanto, si la interpretación de los resultados es tan crucial como la precisión de la predicción, podríamos considerar un balance entre modelos más simples como OLS o modelos de efectos aleatorios y XGBoost.

Además, aunque los modelos SAR y SEM no superaron a XGBoost en precisión, podrían ofrecer insights útiles sobre la estructura espacial de los datos, algo que XGBoost no captura de manera inherente. Por lo tanto, si la dimensión espacial es de interés, valdría la pena explorar estos modelos más a fondo o combinar sus fortalezas con las de XGBoost.

Analisis y Hallazgos

Desarrollar una breve descripción de los 6 – 10 principales hallazgos de: ## a. EDA

1. Tipo de Vehículo: Los montos de reclamación varían según el tipo de vehículo. Los vehículos tipo SUV tienden a tener montos de reclamo más altos, como se muestra en el gráfico de caja.

2. Daños a la Propiedad: Hay una clara distinción en los montos de reclamación cuando hay daños a la propiedad involucrados. Los reclamos con daños a la propiedad tienen generalmente montos más altos.

3. Relación con la Edad: Existe una correlación entre la edad y el monto total del reclamo por accidente. A medida que aumenta la edad, también lo hace la tendencia en el monto del reclamo.

4. Número de Vehículos Involucrados: El número de vehículos involucrados en un reclamo también afecta el monto. Los accidentes con más vehículos involucrados suelen tener montos de reclamo más altos.

5. Correlaciones: Se observa una matriz de correlación que indica la relación entre diferentes variables. Por ejemplo, hay una correlación positiva entre el monto del reclamo del vehículo y el monto total del reclamo.

6. Distribución Geográfica: Los mapas muestran diferencias geográficas en el monto promedio del reclamo y la prima anual promedio de la póliza por estado, así como los límites de cobertura de la póliza por paraguas en diferentes condados dentro de estados específicos.

7. Distribución del Monto Total del Reclamo: La mayoría de los reclamos se concentran en montos más bajos, con una disminución gradual en la frecuencia a medida que aumenta el monto del reclamo, lo que sugiere que los reclamos de alto valor son menos comunes.

b. Modelo seleccionado:

i. ¿Cuáles son las variables que contribuyen a explicar los cambios de la principal

variable de estudio?

Las variables que tienen una contribución significativa al explicar los cambios en la variable de estudio principal, que es el total_claim_amount (monto total de la reclamación), son:

• log_vehicle_claim: El logaritmo del monto reclamado por el vehículo.
• log_injury_claim: El logaritmo del monto reclamado por lesiones personales.
• number_of_vehicles_involved: El número de vehículos involucrados en el incidente.
• property_damage: Si hubo daño a la propiedad.
• police_report_available: Si está disponible un reporte policial.
• incident_hour_of_the_day: La hora del día en que ocurrió el incidente.

ii. ¿Cómo es el impacto de dichas variables explicativas sobre la variable dependiente?

El impacto de estas variables sobre la variable dependiente se observa de la siguiente manera:

• Los montos de reclamo tanto para el vehículo como para lesiones personales, cuando se transforman al logaritmo, tienen una relación positiva y significativa con el monto total de la reclamación.
• El incremento en el número de vehículos involucrados tiende a disminuir el monto total de la reclamación.
• La presencia de daños a la propiedad y la disponibilidad de un reporte policial tienen efectos mixtos y no consistentemente significativos en los diferentes modelos.

iii. ¿Los resultados estimados del modelo seleccionado son similares a los otros modelos

estimados? ¿Cuáles son las diferencias?

Comparando los resultados del modelo seleccionado con otros modelos estimados, encontramos que:

• El modelo OLS (Ordinary Least Squares) proporciona un Multiple R-squared de 0.8588, lo que indica una alta proporción de la varianza explicada por el modelo.
• Los modelos SAR (Spatial Autoregressive) y SEM (Spatial Error Model) incluyen la autocorrelación espacial en los errores, lo que no parece mejorar significativamente el ajuste del modelo en comparación con el OLS estándar, basado en la similitud de los valores AIC.
• Los modelos de aprendizaje automático, como XGBoost, Árbol de Decisión y Bosque Aleatorio (Random Forest), tienen diferentes métricas de rendimiento como RMSE, que varían en magnitud. El modelo XGBoost tiene un RMSE relativamente bajo en comparación con otros modelos, indicando un mejor rendimiento.
• La Red Neuronal muestra una aproximación gráfica de cómo las variables de entrada se relacionan con la salida, pero no se proporciona una métrica de rendimiento específica para la comparación directa.

Los resultados sugieren que los montos de reclamaciones del vehículo y las lesiones personales son factores significativos en la predicción del monto total de la reclamación, mientras que los otros factores tienen un impacto variable y a menudo no significativo. El modelo OLS parece tener un buen ajuste con los datos, aunque los modelos de aprendizaje automático, particularmente XGBoost, podrían ofrecer un rendimiento predictivo más robusto.